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We study the decay of vortices in Bose-Einstein condensates at finite temperatures by 
means of the Zaremba-Nikuni-GrifFin formalism, in which the condensate is modelled by a 
Gross-Pitaevskii equation, which is coupled to a Boltzmann kinetic equation for the thermal 
cloud. At finite temperature, an off-centred vortex in a harmonically trapped pancake- 
shaped condensate decays by spiralling out towards the edge of the condensate. This decay, 
which depends heavily on temperature and atomic collisions, agrees with that predicted by 
the Hall-Vinen phenomenological model of friction force, which is used to describe quantised 
vorticity in superfluid systems. Our result thus clarifies the microscopic origin of the friction 
and provides an ab initio determination of its value. 



INTRODUCTION 



The dynamics of Bose-Einstein condensates at finite temperature presents an interesting prob- 
lem in the study of ultra cold Bose gases. In most experiments, such systems are only partially 
condensed, with the non-condensed thermal cloud providing a source of dissipation and leading to 
damping of structures, such as collective modes [l|, 0, 0, [J, solitons [H] and vortices 0]- Sev- 
eral approaches have been developed to describe these systems, including generalised mean field 
treatments 0, [13; lI^HjIlS]' number-conserving approaches 14, 15, 16], classical field theory 
12, 18, 19], and stochastic approaches 2^, 21, 22], as recently reviewed by two of the authors 23], 
who give a more complete list of references. Although the underlying theory is well understood, 
the implementation of models which can be actually solved in specific contexts has proven to be 
a considerable challenge, with the majority of treatments to date assuming the thermal cloud is 
homogeneous and static. 

In this paper we use the formalism of Zaremba, Nikuni and Griffin (ZNG) [l3]- The ZNG theory 
is a kinetic approach in which a generalised Gross-Pitaevskii (GP) equation for the condensate 
order parameter is coupled to a Boltzmann equation for the thermal cloud. These equations have 
already been solved numerically 24], and applied to the study of collective modes 25, 2^, 27], the 
hydrodynamic regime 28, 22] and the decay of dark solitons [soj, demonstrating good agreement 
with experiments. 

The primary aim of this paper is to apply the ZNG model to quantised vortices at finite 
temperatures. In a harmonically-trapped condensate at zero temperature a vortex will precess 
within the condensate, following a trajectory of constant radius, for which the energy remains 
constant. At finite temperature the presence of dissipation leads to the vortex minimising its 
energy by moving towards the surface, where it eventually leaves the condensate and disappears. 
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This decay has been observed experimentally 0, 0], but the precise theoretical modelling of this 
process has proven rather challenging. 

Our approach should be contrasted with earlier work. Dissipation based on the scattering 
of single-particle excitations from the mean field potential at the vortex core was considered by 
Fedichev and Shlyapnikov 31 1 and was subsequently extended to treat vortex lattices 32, 3^. The 



problem of vortex nucleation was then considered by Penckwitt et al. 37[ using an approximation 



in which the thermal cloud is treated as static. Schmidt et al. 3j] were the first to apply the 
classical field method to the problem of vortex decay, while more recent applications of the method 
addressed vortex dynamics in quasi- two-dimensional systems IE, 36]. Duine et al. [s^ considered 
the dynamics of a straight vortex line using a closely-related method based on a stochastic Gross- 
Pitaevskii equation [2l| and derived a stochastic equation of motion for the position of the vortex 
core; a similar equation was also used by Sasik et al. to numerically simulate the motion of 
an isolated vortex in a uniform box with periodic boundary conditions [39]. Compared to the 
literature cited above, our work represents the first microscopic simulations which fully account for 
the dynamics of an inhomogeneous thermal cloud. 

Our work, however, also has a second motivation. Although the topic of quantised vorticity in 
superfluids is interesting per se (as shown by the number of recent vortex experiments in atomic 
Bose-Einstein condensates), it also has implications in the subject of quantum turbulence [40*]. 
Current work on turbulent superfluid ^He and ^He-B, for example, is concerned with the extent to 
which turbulence in these systems differs from that found in ordinary classical fluids [4ll . [13, SB 
45l. l46l. 1471. l48l. 1491] . What makes quantum fluids attractive from the point of view of understanding 
the principles of turbulence, is the existence of various forms of dissipation which are distinct 
from ordinary viscosity. At sufficiently low temperatures, kinetic energy can be dissipated into 
sound waves, that is phonons [50(], via a Kelvin wave cascade 5]| or via vortex reconnections 
[s^ ]. At higher temperatures, the friction force between the superfluid and the normal fluid 
component can change the nature of the turbulent kinetic energy cascade. For example, in the 
classical turbulence scenario [53], the Richardson-Kolmogorov inertial cascade is limited at large 
wavenumbers where viscous dissipation destroys the small scales, whereas in superfluid ^He-B at 
relatively high temperatures the friction can limit the inertial cascade at small wavenumbers 55]. 

The mutual friction between a superfluid and a normal fluid is one of the most intricate is- 
sues of superfluidity jsB]; in particular, the existence of a transverse force on quantised vortices 
parametrised by a dimensionless temperature dependence quantity called a' has been controversial 
57 . 58 . 59 . 6o| . Studying vortex motion in an atomic BEC at flnite temperatures, and interpreting 



the results from the point of view of vortex dynamics, allows one to compute the friction force 
directly from flrst principles. This work on an atomic quantum fluid thus provides insights into 
this important problem which cannot be obtained as readily with superfluid helium. 

This paper is structured as follows: Sec. HI] briefly reviews the ZNG theoretical model and 
numerical implementation. Sec. IIIII presents our main findings on vortex decay (Sec. IIII Ap . high- 




FIG. 1: (Colour online) Three-dimensional density isosurface showing an off-centred vortex in a pancake- 
shaped condensate. 
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lighting the dependence of the friction coefficients on system parameters (Sec. IIIIBp . the role of 
collisions between the atoms (Sec. IIIl"Cl) . and the effect of thermal cloud rotation (Sec. IIIID]) . Sec. 
IIVI presents some concluding remarks and briefly discusses the consequences of our analysis on the 
motion of vortex lattices. 



II. THEORY 
A. ZNG formalism 



Following Refs [10| and [2j], the second-quantised field operator ip{r, t) can be split into conden- 
sate and thermal contributions. Making use of Bose broken symmetry, the full operator is written 
as ip{r,t) = ^(r,i) +i/)(r,f), where ^{r,t) = {ip{r,t)) is the condensate wavefunction (angular 
brackets denote an expectation value) and il){r, t) is the noncondensate field operator. Starting 
with the Heisenberg equation of motion for ip(r,t), one eventually arrives at the following pair of 
equations: 

5^- / n^v^ 



ih—=^-^ + V + gnc + 2gn-iRj^, (1) 

^+P.Vf-VU-Vpf = Cu + C22. (2) 

ot m 

Eq. ([1]) is a generalised Gross-Pitaevskii (GP) equation for the condensate wavefunction ^(r,t), 
and has been obtained in the so-called 'Hartree-Fock-Popov' approximation whereby the static 
value of the 'anomalous' average, which is responsible for certain many- body effects (eil,^,^,!^!]- 



is ignored. Eq. ([2]) is a Boltzmann equation for the thermal cloud phase space density f{p,r,t), 
with thermal energies calculated in the Hartree-Fock approximation 0. The condensate density 
is defined as ric = while the thermal cloud density is obtained from / by means of the 

momentum integral n = f dp/h^ f. The mean field interactions between atoms is parameterised by 
g = Anh'^a/'m, where m is the atomic mass and a is the s-wave scattering length. The thermal atoms 
experience an effective potential given by U{r) = V{r)+2g{nc-'rh), where V{r) = m(w^r^+w^2:^) /2 
represents the external trap. 

The terms involving gric and gh in Eq. ([T|) and VC/ in Eq. ([2]) represent mean field coupling 
between atoms in the condensate and the thermal cloud. This coupling is a source of dissipation 



for the system, and gives rise, for example, to Landau damping of collective modes j66l. 1671. l68l. 1691]. 
The ZNG model also includes 'coUisional integrals' C22 and C12, which respectively denote binary 
collisions between noncondensate atoms, and between condensate and noncondensate atoms. They 
are given by: 



C22 = I dp2dp^dp,6{p+P2-p3-p4)S{e+€2-€s-e,)[{l+f){l+f2)f3h-ff2{l+f3){l+h)], 

(3) 

2g'^nc f 

C12 = ^2TTyh^ J ^P2dp3C^P4<5(mi'c +P2 -P3 -P4)<5(ec + £2 - £3 - e4)['5(p -P2) - '5(P -Ps) - '5(P -P4)] 

x[(l + /2)/3/4-/2(l + /3)(l + /4)], (4) 



where / = f{p,r,t) and fi = f{pi,r,t). In the above expressions, delta functions enforce 
momentum and energy conservation in the collisions, where e = /{2m) + U \s the thermal 
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atom energy (in the Hartree-Fock limit), Ec = mv'^/2 + is the local condensate energy, with 
Vc = h{'^*V'^ — ^'V^*)/(2imnc), and /ic is the chemical potential. 

The Ci2 term (jH) involves those collisions between condensate and thermal atoms which lead to 
a transfer of atoms between condensate and thermal cloud. This term is thus related to the source 
term —iR^ appearing in Eq. ([T|) via 

If R is positive (negative), there is a net local flux of atoms out of (into) the condensate. 



Numerical methods 



The methods used for our numerical simulations are discussed more fully in Ref. [2J], and are 
only briefly reviewed here. 

Firstly, we must generate a suitable initial state for the simulations which consists of a con- 
densate in equilibrium with a thermal cloud at temperature T. The condensate wavefunction can 
be obtained by an imaginary time propagation (t — > —it) of ([T]) with R = Q. The thermal cloud 
density n = is first set to zero, and an approximate thermal cloud potential [/ ~ V + 2gnc is 
constructed from the self-consistently determined condensate wavefunction. This yields the ini- 
tial thermal cloud density nQ{r) = (1/A^) (73/2(^)1 where A = (Inh^ /mkBTy^'^ is the thermal de 
Broglie wavelength, z{r) = exp{/3[//c — U{r)]} is the local fugacity and fic is the condensate chemi- 
cal potential. This density is then used in ([T]) to obtain an improved condensate wavefunction and 
the procedure is iterated until a self-consistent solution of both the condensate and thermal cloud 
densities is obtained. (For more details see [j^].) 

The vortex state of interest in our simulations can be obtained from this equilibrium state by 
multiplying the condensate wavefunction by the phase factor exp[i5(r)], where S{r) = arctan((7/ — 
yo)/{x — xq)) is the phase profile associated with a straight vortex located at (2:0,2/0) • The GP 
equation is then solved again in imaginary time for a short period, until a vortex is fully formed 
in the condensate and most short time scale transients in the initial configuration (for example, 
phonon excitations) have damped out. 

The state generated in this way is a quasi-equilibrium state containing one vortex whose sub- 
sequent dynamical evolution is of interest. The generalised Gross-Pitaevskii equation for the con- 
densate wavefunction ([T]) can be solved readily in real time using standard methods [t^], in our 
case a split-operator Fast Fourier Transform approach on a 3D Cartesian grid. The thermal cloud 
is described by a swarm of classical test particles, which move in an external potential U{r, t), and 
which, during a timestep, can collide with each other or with the condensate. The probabilities 
of particle collisions are chosen so that they correspond to a Monte Carlo evaluation of the col- 
lision integrals in ^ and (H]). Together with the Newtonian dynamics of the test particles, this 
procedure is equivalent to solving the collisional Boltzmann equation ([2]). The C12 probabilities 
are then summed according to ([5]) to determine R{r,t), which appears in the GP equation ([1]). 
Since R{r, t) is a non-Hermitian term, it leads to change in the normalisation of the wavefunction, 
corresponding to condensate growth or loss, which is accompanied by a compensating removal or 
creation of thermal particles. An essential ingredient in the simulations is the evaluation of the 
thermal cloud density h{r,t), which appears in both ([T]) and ([2]) (through the effective potential 
U) . This is achieved by appropriately binning the thermal particles and then convolving the binned 
distribution with a Gaussian in order to obtain a smoothly varying potential. 
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FIG. 2: (a) Position of the vortex (in units of a± = y^h/muj±), as it spirals out of the condensate 

at T = 0.7Tc. (b) Radial position of the vortex r„ — \/ x'^ + y'^ as function of time (in units of for 
T = O.STc (solid black line), T = O.GTc (dashed black), and T ^ 0.7Tc (solid grey). 



III. RESULTS 



The first simulation to be described is for a pancake-shaped condensate with A'^ = 10^ ^^Rb 
atoms with trap frequencies u!± = 2it x 129 Hz and uJz = ^/8u!±. For these parameters the (ideal 
gas) critical temperature is Tc = 177 nK. This geometry has the advantage that the radius of the 
condensate in the axial direction is much smaller than in the radial, so that the vortex remains 
relatively straight throughout its motion. This simplifies the analysis considerably, as the vortex 
dynamics can be characterised by just its x and y coordinates. 



A. Vortex decay 



Following the initial preparation phase, the vortex is centred on x^(0) ~ 1.3a±, y^(0) = 0, as 
illustrated in Fig. ([1]) , where a± = \J h/muj±_ is the harmonic oscillator length in the radial direction. 
The subsequent vortex position r^(t) = {xy{t),yv{t)) is tracked by finding the corresponding local 
density minimum in the z = plane by means of a quadratic interpolation between grid points. 
This procedure breaks down when the vortex leaves the bulk of the condensate into the very low 
densities beyond the edge, so no results are shown when this happens. This condition can be taken 
as the point at which the vortex "disappears" from the condensate, and would correspond in the 
experimental context to the density contrast of the vortex being below the detection limit. Since 
the initialisation process induces a centre-of-mass motion of the condensate, the following analysis 
depicts the vortex position relative to the centre-of-mass position. 

Simulations of the GP equation for T = reveal that the vortex precesses in a circular path 
around the condensate, following a trajectory of constant energy as would be expected for a non- 
dissipative system. This well-known precessional behaviour can be understood as arising from the 
non-uniform density of the condensate, which means that the energy of the vortex is a function of 
its radial position. Hence there is an effective Magnus force, proportional to the gradient of the 
energy and directed radially, with in turn induces the azimuthal vortex motion 71 1. This can be 



described quantitatively by means of a time-dependent variational method 721. 173l.l74i|. In addition 



to the motion described above, the acceleration experienced by the vortex in its circular trajectory 
can in principle lead to the emission of sound waves. However, for the harmonic confinement 
being consider, any emission of sound waves is followed by reabsorptioon, and no net dissipation 
occurs. There is nevertheless some modulation of the vortex trajectory arising from the dynamical 



interaction between of the vortex and sound waves 751 ]. 



In contrast to this dissipationless motion, the vortex can lose energy at finite temperatures 
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due to interactions with the thermal component, and as a result, it moves radially towards the 
condensate edge. The resulting spiral trajectory is illustrated in Fig. [2ja) for a temperature of 
T = 0.7 Tc. Fig. [2]^b) shows the outward relaxation of the radial position r„ = |r^| = yxj+yj 
as a function of time. The three curves represent different temperatures, and demonstrate, as one 
would expect, that the relaxation rate increases with temperature. 

As well as monitoring the vortex position, it is instructive to study the evolution of the densities 
during the simulation. Fig. [3] shows the condensate (top images) and thermal cloud densities 
(bottom) (cross-sections at z = 0) for various times, and for a temperature of T = 0.7 Tc. The 
density dip at the vortex core (dark blue) is evident in the condensate images (top), and executes 
the spiral motion of FiglJJa). The bottom images show the thermal cloud density; as this is much 
smaller than the condensate density, a different colour scale is used here for clarity. The thermal 
cloud density is largest near minima of the effective potential U{r,t), leading in the present situation 
to the circular ring of peak thermal cloud density (red) at the edge of the condensate where ric — > 0. 
However, the thermal cloud also tends to "fill in" the vortex core since the lower condensate density 
at this position gives rise to a dip in the effective potential. The resulting peak in the thermal 
cloud density tends to follow the core as it precesses and spirals out. We note that the stabilisation 
of the vortex by the thermal cloud suggested in 76[] is not observed. It arises when the thermal 
cloud is treated as static [t^, but does not occur when the dynamics of the thermal cloud is taken 
into account. 

The radial position of a vortex close to the trap centre exhibits a near-exponential growth in 
time. This is evident from the inset of Fig. HI which plots the time evolution of on a logarithmic- 
linear scale for T = 0.5 Tc. Deviations from this simple exponential behaviour are observed at 
later times when the vortex approaches the edge of the condensate; similar behaviour is found for 
other temperatures. To quantify the relaxation we thus fit r^(i) over < t < 50 to the function 
r^{t) = ro exp(7t). The resulting values of 7 are plotted in the main panel of Fig. H] with the black 
circles. 

The departure from a pure exponential behaviour can be interpreted as a position-dependent 
relaxation rate 7 which is enhanced by the local maximum in the thermal atom density near the 
edge of the condensate (see Fig. [3]). To check this interpretation, we performed simulations in 
which the vortex starts nearer to the centre, e.g. at tq — 0.65a_L, instead of the initial position 
ro — 1.3aj_ used earlier. The resulting values of 7 are plotted as the open circles in Fig. HI showing 
smaller values consistent with the lower thermal cloud density closer to the centre of the trap. 
However, the difference is small, showing that the exponential decay approximation is good until 
quite near to the condensate edge, where the larger thermal cloud density leads to a more rapid 
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FIG. 3: (Colour online) Density cross-sections of the condensate (top row) and thermal cloud (bottom row) 
at z = for T = 0.7 Tc, and at the times (a) uj±t ~ 6, (b) 12 (c) 18, and (d) 24. The colours range from 
brown/red (high density) to dark blue (low density), with different scales for the condensate and thermal 
cloud densities. 
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FIG. 4: (Inset) Lin-log plot of the vortex radial position as a function of time, for T = 0.5 Tc (solid line). 
The dashed line is an exponential fit, r^{t) = tqC^*, to the data over Q <t < 50. The main figure plots the 
resulting values of 7 for for an initial vortex position of rp ~ 1.3 (solid circles) and ro — 0.65 (open circles). 
For comparison, the solid and dashed lines plot the results of FS [sij and DLS [s^ respectively. 



relaxation (inset of Fig. S]) . 

It is of interest to c omp are our results to existing analytical pre dictions obtained by Fedichev 
and Shlyapnikov (FS) [3l|, and Duine, Leurs and Stoof (DLS) [381]. Both FS and DLS actually 
quote the timescale t for the decay of a vortex from position Tmin to 

'"max) but this can be converted 

to a decay rate using r = ln(rmax/'^min)- In both of these works the decay of the vortex is 
found to be exponential, with a rate 7 oc T but with different proportionality coefficients due to 
the different approximations made in the two theories (see below for the role of the noise term in 
the DLS analysis). These rates are displayed in Fig. Hlby solid (FS) and dashed lines (DLS). It is 
apparent that the rates found by FS are comparable to ours at the lower temperatures, but differ 
significantly at higher T due to the stronger (approximately quadratic) temperature dependence 
found in our simulations. 

FS assumed a uniform condensate in a cylindrical container, and modelled the decay solely as 
the result of mean-field interactions. The DLS study on the other hand, includes the important 
C12 collisional coupling between the condensate and thermal cloud. With the aim of obtaining 
analytical results, they approximate the condensate density profile by a Gaussian which is rea- 



sonably accurate for the most relevant region near the center of the trap 77|]. In this regard, it 
should be noted that such an approximation has been shown to produce correct results for the 
frequencies of collective modes even for Thomas- Fermi condensates 78|. Our present simulations, 
which include both mean field and collisional coupling mechanisms, enable us to assess their relative 
importance, which will be discussed in Sec. IIII CI More importantly, our simulations differ from 
these approaches in that they are actually performed for a dynamical thermal cloud, with both the 
condensate and thermal cloud densities determined self-consistently during the simulations. 

For completeness, we should however make two additional remarks regarding the DLS approach. 
Firstly, in their preceding work [t^, Duine and Stoof argued that enhanced damping of collective 
modes at higher temperatures could be related to the position dependence of the self-energy (and 
hence of the damping term —iR^), which was ignored in their analytical treatment based on a 
volume average over the size of the condensate. In the present context, this could partly account 
[7^ for a deviation of the computed damping rate ^/oj± at higher temperatures from the linear 
behaviour seen in Fig. HI Although such simulations have not been performed to date, the resulting 
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corrections are likely to be smaller than the observed disagreement, whose origin we believe lies 
primarily in the dynamics of the thermal cloud. 

On the other hand, the DLS analysis is actually more general than ours in that it contains an 
additional noise term in the equation of motion of the vortex. This provides stochastic "kicks" to the 
vortex and the ensuing Brownian motion allows for the migration of a vortex away from the centre 
of the trap. Were it not for the noise, a centred vortex would have an infinite lifetime. This would 
be the case in an exact application of the ZNG theory. However, its numerical implementation in 
terms of discrete test particles does introduce statistical fluctuations in the thermal cloud density 
which plays the role of noise. Thus we indeed find in simulations of a centred vortex a finite, 
albeit long, lifetime (wxr ~ 1000 for the relatively high temperature of T = 0.7 Tc as compared to 
uj±T ~ 120 for r^(0) ~ 1.3a±). This long lifetime, however, should not be taken seriously since it 
depends on the actual number of test particles used in the simulations. The simulation nevertheless 
makes clear that this 'numerical noise' is of secondary importance at larger radii where the direct 
coupling to the thermal cloud is the dominant dissipative effect. Although it is something to be 
checked, it is unlikely that the stochastic term makes a significant contribution to the spiralling 
out of a vortex when it is located far from the trap centre. 

Finally, it is worth remarking that our computed decay rates (from approximately 0.5 s^^ to 
3 in the range T/Tc = 0.4 to 0.6) are in order-of-magnitude agreement with the decay rates 
observed for a vortex lattice [3] (approximately 0.3 s~^ to 3 s~^ over the same relative range). 



B. Friction coefficients 



In order to further understand the origin of this exponential decay, it is instructive to consider 
the two-fluid hydrodynamics model used to describe superfluid liquid helium [io|. In this context, 
dissipation arises from the interaction between the quantised vortices and the thermal excitations 
(phonons and rotons) which form the normal fluid [s^ . 81]. Since the radius of the superfluid 



vortex core is much smaller than the typical separation between vortices or any other length scale 
of interest in the flow, the vortex is described in parametric form as a three-dimensional space 
curve s = s(^,t), where ^ is the arclength. The resulting equation of motion [s^] is: 



ds 
dt 



+ Vi + as' X (v„ - - Vj) - a's' x [s' x (v„ - - Vj)], 



(6) 



where s' = ds/d^ is the unit tangent along the vortex at the position s, is any imposed superfluid 
velocity, v„ is the normal fluid velocity, and Vj is the self-induced velocity of the vortex arising 
from its own curvature, the presence of other vortices, and any inhomogeneity of the fluid. The 
first two terms on the right hand side of Eq. ([6]) state that at T = the vortex is advected by the 
local superflow + Vj. The remaining two terms reflect the fact that, at nonzero T, the normal 
fluid streaming past the vortex core exerts a force per unit length whose intensity is controlled by 
the temperature-dependent friction coefficients a and a' [53l . 

In turbulent helium, the friction coefficients a and a' play a key role because they account 
for the mutual friction between the superfluid and the normal fluid and control the transfer of 
energy between the two fluids at various length scales and hence the nature of the inertial range 
cascade [i^, Isp]. The first coefficient, a, describes dissipative effects and leads, for example, to 
the shrinking [53] of a vortex ring or the damping of a Kelvin wave 85] in a normal fluid at rest. 
In the case of rotating helium, a determines the attenuation of second sound waves, so it allows 
the experimentalist to determine the density of vortex lines ^8^. The second coefficient, a', is not 
dissipative, and, in rotating helium, splits a second sound resonance (besides the classical rotational 
splitting) in a suitably designed cavity 87]. It has been argued 89, [o^] that the ratio of inertial 
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FIG. 5: (a) Precession angular frequency of the vortex, as a multiple of the radial trap frequency uj^. (b) 
Friction coefficient a as a function of relative temperature T jTc, for TV — 10"* and a pancake trap geometry 



and dissipative forces, which is called the Reynolds number in the case of ordinary turbulence 54i |. 
is simply a/(l — a') in the case of quantum turbulence. 

In our case of a pancake-shaped condensate the vortex line remains approximately straight, so 
we can then replace s by the vortex position r^,, while s' = z. Both the condensate and thermal 
cloud are stationary, and so = v„ = 0. This just leaves Vj, which in this case corresponds to 
the azimuthal vortex motion induced by the inhomogeneity of the condensate at zero T (hence, 
Vj = Vi(j)). Rewriting Eq. ([6]) in cylindrical polar coordinates gives for the azimuthal component 



(1 - a ) — , 



(7) 



where lo^ = (j)^ is the vortex precession frequency at finite T and the dot denotes a time derivative. 

To a first approximation we ignore the mutual friction coefficient a'. We then find oj^ = Vi/r^, 
and thus obtain for the radial component 



dry 
dt 



(8) 



If a and u)^ are constant, then one simply recovers the exponential behaviour discussed in the 
previous section, with 7 = aujy. 

In actual fact, a and uJi, would not be expected to be constant during the course of the vortex 
decay. The condensate and thermal cloud are both non-uniform (implying a dependence of a on 
position), while oj^ is only constant near the centre, and increases as the vortex approaches the 
edge [71|. However, from the evidence of Fig. H] (inset), Eq. ([8]) is a good approximation for times 
when the vortex is close to the centre where the densities are relatively uniform. 

Values of a can therefore be estimated from 7 in Fig. HI The other necessary ingredient is 
bjy, which can be calculated as a function of time using uj^ = Xyi/y — i/vXy, where numerically 
the time derivatives are calculated using central differences. Due to errors in locating the vortex 
position, there are large fluctuations in the calculated tUy. To obtain the corresponding values we 
thus average results over 2 — 3 orbits of the vortex, corresponding to times < uj±t < 50. The 
resulting mean values are plotted in Fig. [5] (a). These results are also used for calculating the 
values of a presented in Fig. [5] (b). The observed increase in a with rising T is similar to the 
behaviour observed in liquid ^He 84 1 and ^He 53, [s^], and is consistent with results obtained 
using the related but simpler phenomenologically damped GP equation[91]. 

Let us now investigate the role of the a' coefficient in Eq. ([6]). Experiments in liquid helium 
show a small effect, and there has been some controversy in the literature about this transverse 
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component of the friction force 57, 58, 5^, 60]. This issue can be addressed within our simulations 
by comparing the "dynamic" thermal cloud frequency iv^ to a "static" value Wgt found using a 
GP simulation in which the dynamics of the thermal cloud is ignored, that is, it retains its initial 
equilibrium form for all times. In this static thermal cloud approximation, the thermal cloud exerts 
a time-independent mean-field potential on the condensate, and its effect can be identified with 
the Vi/r^ term in Eq. ([7]), giving the simple relation for a' = (1 —LOy/uJst)- We find that w^, and cogt 
are equal to within 2 — 3 %, showing that the changes observed in ujy in Fig. [5] are almost entirely 
due to the effects of the static thermal cloud potential on the condensate density profile, and not 
to "real" dynamical effects of finite temperatures. The errors in measuring the vortex precession 
frequency are such that we cannot confidently extract a value for a' , although our simulations 
indicate that \a'\ < 0.02 throughout the measured temperature range < T/Tc < 0.8. This agrees 
with recent results of Berloff and Youd 92] obtained by means of classical field theory. 

We also explore the dependence of these parameters on the total number of atoms, N. Fig. 
[6] shows results for between 10^ and 10^, and a fixed value of T = 0.5 Tc. The decay rate 7, 
plotted in Fig. [6] (a), tends to decrease with increasing N. The vortex precession frequency shown 
in Fig. [6] (b) also has this decreasing trend, in agreement with what one would expect from GP 
solutions. These decreasing trends approximately cancel when calculating a = ^/tOv, leading to 
variations having no clear dependence on A^. The scatter reflects the uncertainty in the calculation 
of a and suggests that this parameter is approximately independent of atom number. 



X 10" 




X 10 



FIG. 6: (a) Decay rate 7, (b) precession frequency oj^, and (c) a, as a function of total number of atoms N, 
for temperature T/Tc = 0.5 and a pancake geometry. 



C. Collisionless simulations 

In general, damping in the ZNG formalism arises from the coupling of the condensate to the 
thermal cloud by means of mean field interactions and C12 collisions. In order to explore the 
relative importance of these two contributions, we have performed simulations where collisions are 
not present, so C12 = C22 = in ([2]), and the only source of dissipation is mean- field coupling. 
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FIG. 7: Radial position of the vortex vs. time for T — 0.6 Tc, and coUisional (black line), and coUisionless 
(grey line) simulations. 



Physically, this dissipation is a form of Landau damping whereby the motion of the vortex core 
through the thermal cloud generates thermal excitations [6^ . In Fig. [7] the radial position for the 
coUisionless simulation at T = 0.6 Tc is shown as the grey line, and compared to the collisional 
result in black. Without collisions one obtains a much slower decay, highlighting the crucial role 
of collisional damping in the simulations, a conclusion which has been numerically verified over a 
broad range of temperatures. The near linear variation of r„ with time seen in Fig. [7] indicates that 
an exponential function would be a poor fit to the data in this case. In view of Eq. it would 
appear that the mean- field contribution to a must be a decreasing function of but we have no 
simple explanation for this. 



D. Rotating thermal clouds 



We also consider the case where, instead of being stationary initially, the thermal cloud under- 
goes solid-body rotation around the z-axis with angular frequency f^th- The thermal cloud velocity 
at the vortex core is then Vn = fith'"i;0- Hence r^th > represents rotation in the same sense as 
the vortex precession, while f^th < indicates rotation in the opposite sense. Using Eq. ^ then 
yields: 



a;^, = (1 - a') — + a'r^th, 



and. 



r„ = ro e 



(9) 



(10) 



For these simulations we start with the equilibrium condensate and thermal cloud distributions 
evaluated at T = 0.7 Tc and Jlth = 0. A rigid body rotation of the thermal cloud is then imposed 
by adding v„ = Q^i^r^cf) to each atom's velocity. It should be noted that the thermal cloud will now 
no longer be in "equilibrium" since there is a centrifugal effect which tends to expand the cloud. 
This initial outward expansion leads to an oscillation in the radial direction and, since angular 
momentum is conserved, a corresponding oscillation of the cloud's angular velocity. 

The vortex decay curves are plotted in Fig. ([8]) for different Qth, and show that the decay rate 
increases for rotations opposite to the vortex precession direction (Q^h = —0.2), but decreases when 
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FIG. 8: (Colour online). Time evolution of vortex radial position for T — 0.7 Tc and a rotating thermal 
cloud. The different curves represent varying thermal cloud rotation rates, with f^th = —0.2 (black, top), 
f^th = (red), Slth = 0.2 (green), and f^th = 0.37 (blue, bottom). 

they rotate in the same direction (f^th = 0.2 and Oth = 0.37). This is consistent with the expected 
behaviour from Eq. (|1U|). To study the problem more quantitatively, we again fit exponentials of 
the form ae^^ to the decay curves over < u!±t < 50. The values of b for the different 0th ^ire 
plotted in Fig. [H The straight line is the expected result from (fTOl) . b = a{uJv — ^th), where a 
and ojy are taken from the results of the f^th = simulations found earlier. Our results are in 
quite good agreement with this behaviour, although some small discrepancies are apparent. These 
could be due to the oscillations in Q^h noted earlier, whereas Eq. ([10]) assumes that Jlth is strictly 
time-independent throughout the precessional motion of the vortex. 

0.015 



0.01 
0.005 



-0.2 n 0.2 

"th 

FIG. 9: Value of decay rate 7 for different rotation rates of the thermal cloud, ilth- The solid line plots 
the function 7 = 7o(l — ^th/i^vo) (where 70 and lOvq are the decay rate and the precession angular velocity 
respectively for a non- rotating thermal cloud), which is the expected dependence from Eq. (|10|1 . 
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IV. CONCLUSIONS 

In summary, we have studied the finite temperature dynamics of a single vortex in a partially- 
condensed ultra-cold Bose gas. Our methodology and detailed simulations provide several advan- 
tages over previous studies. Firstly, our simulations include the full effects of the trapping potential 
through the self-consistent determination of condensate and thermal cloud in the initial state. This 
results in a more realistic model as compared to those using the approximation of uniform den- 
sities 3l|, l37l . ISSi ] . The inclusion of non- uniform densities accounts more realistically for both the 
dynamics of the vortex and the positional dependence of the dissipation. Secondly, our model 
includes both mean-field j3l| and collisional [ss] damping, and allows us to compare the relative 
importance of the two mechanisms. Thirdly, the thermal cloud is not assumed to be static as in 
earlier treatments, but is treated dynamically on the same footing as the condensate. This more 
refined treatment negates suggestions that the thermal cloud can act as a "pinning potential" , sta- 



bilising the vortex 76|]. In addition, it has allowed us to observe the change in damping when the 
thermal cloud is moving relative to the condensate as, for example, when undergoing a rotation. 

We also compared our results for the vortex relaxation rate to those of other studies and found 
some significant differences, particularly with regard to its temperature dependence. Furthermore, 
by comparing the trajectory of vortices with the predictions of phenomenological vortex dynamics 
equations, we were able determined the mutual friction coefficients from first principles. 

Our approach can also be extended to study the role of a dynamical thermal cloud on vortex 



lattice dynamics |93l . 194 l95l | , thereby complementing and extending existing work [32, l33|, l37|, l7a, 
96, 97, 98, [9^. To illustrate this possibility, we conclude by briefly reporting on some preliminary 
results for the decay of vortex lattices. We consider the evolution of two different vortex lattice 
configurations shown in Fig. [TOj Both vortex arrays initially contain seven vortices (left images), 
however they differ in the way the vortices are arranged. The first array (top images) consists of 
one vortex at the centre of the condensate and a ring of six vortices around it, whereas the second 
array (bottom) consists of a ring of seven vortices with no central vortex. These arrays rotate 
in the laboratory frame, and at finite temperatures, the effects of dissipation lead to the gradual 
disappearance of the off-centred vortices, one by one. For simulations performed at T = 0.7Tc, we 
find that the array with all vortices arranged in a ring decays faster: after the first six vortices 
have decayed, the system is left with a single off-centre vortex which moves relatively rapidly to 
the edge and disappears. This whole evolution occurs on a time scale u!±t ^ 150. The decay rate is 
in order-of-magnitude agreement with measurements performed with a bigger lattice 0] after the 
latter are extrapolated to our value of T/Tc- 




FIG. 10: (Colour online). Vortex lattice decay at T = 0.7 Tc- The top row shows the time evolution (left to 
right) of a lattice with a vortex initially at the centre; the bottom row shows the evolution of a lattice with 
the same number of vortices but having a ring configuration. 

In contrast to this, the lattice with the central vortex reaches a point where a single metastable 
central vortex remains after the other six have been shed. This vortex also eventually decays, 
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but our simulations suggest that the decay occurs on a much longer timescale. Some numerical 
experiments we have performed for configurations with no initial central vortex have exhibited a 
similar metastable behaviour. If, during the initial part of the evolution (in which the vortices 
move irregularly), a vortex ends up sufficiently close to the centre, it can become "stuck" nea r the 
centre while the other vortices are shed. These observations are in agreement with reports lod l 
that the decay time of the last vortex is much longer than that of the initial vortex array. 

We stress, however, that a more accurate treatment of the evolution of a metastable central 



vortex requires the explicit inclusion of stochastic noise to provide a "kick", as discussed in 38l. l39l|. 
Such a term was included in a recent discussion [o^, however the thermal cloud was still treated 
as static. A combination of the stochastic Gross-Pitaevskii equation and the quantum Boltzmann 
equation may thus be needed to provide a more complete description of this particular situation. 
However, we expect that most other cases can be modelled extremely well by the Zaremba - Nikuni- 
Griffin approach. In particular, it would of interest to see if more detailed calculations of vortex 
lattice decay would be consistent with experimental observations [T]. Other applications might 
include the study of vortex lattice excitations (Tkachenko modes) [lOll ]. and the dynamics of bent 
vortices in elongated condensates fot]. 

This work was funded by EPSRC grant EP/D040892/1 (BJ, NPP, CFB) and by NSERC of 
Canada (EZ). 
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